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We use the dynamical density matrix renormalization group technique to show that the high- 
energy part of the spectrum of a S = 1 Heisenberg chain, placed in a strong external magnetic 
field H exceeding the Haldane gap A, contains edge singularities, similar to those known to exist 
in the low-energy spectral response. It is demonstrated that in the frequency range ui > A the 
longitudinal (with respect to the applied field) dynamical structure factor is dominated by the 
power-law singularity 5 1 " (q = n, u) oc (to — <jJo)~°" ■ We study the behavior of the high-energy edge 
exponent a' and the edge uio as functions of the magnetic field. The existence of edge singularities 
at high energies is directly related to the Tomonaga-Luttinger liquid character of the ground state 
at H > A and is expected to be a general feature of one-dimensional gapped spin systems in high 
magnetic fields. 
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I. INTRODUCTION 

Studying spectral response is a valuable tool to reveal 
the properties of the strongly correlated ground state in 
interacting electronic systems. One of the paradigmatic 
concepts in physics of one-dimensional (Id) systems is 
the so-called Tomonaga-Luttinger liquid (TLL), 1 which 
for Id systems plays a similar role as the Fermi liquid 
theory for higher dimensions. One of the prominent fea- 
tures of the TLL is the absence of a quasiparticle peak in 
the spectral function; instead, there is a power-law sin- 
gularity, with a non-universal exponent that depends on 
the interaction strength^ 

Apart from Id conductors (e.g., such as carbon 
nanotubes 4 -^) and edge states in fractional quantum Hall 
systems^" 7 - the TLL ground state is expected to exist in 
several other Id systems, particularly in spin chains and 
ladders. For an antiferromagnetic (AF) S = | Heisen- 
berg chain, the TLL character of the ground state can be 
rigorously established from the Bethe ansatz solution, 8 
and it is believed 9 -^ that this picture is valid also for the 
other gapless AF chains with half-odd-integer S. 

For AF Heisenberg chains with integer S and for spin 
ladders, which have a finite excitation gap A, it has 
been argued by analytica l 11 ! 12 and numerica l 13 ! 14 meth- 
ods that in an external magnetic field H > A, strong 
enough to close the spin gap and cause a finite magneti- 
zation, the ground state is also of the TLL type. 

Recently, in a search for Luttinger liquid signatures 
in spin systems, several experimental studies have been 
undertaken i 15 : 16 ' 17 : 18 ! 19 Most of the experimental evi- 
dence is, however, indirect, based, e.g., on the analy- 
sis of the temperature dependence of specific hea t 17 ! 19 
or NMR relaxation rate? 1 ^ Direct detection of the low- 
energy excitation continuum in inelastic neutron scatter- 
ing experiment s] 15 ! 16 is a difficult task, not to mention 



extracting the dynamical exponents from the low-energy 
spectrum. High-energy modes might be easier to study, 
especially with techniques such as electron spin resonance 
(ESR), 2 ^ 21 . 

At the same time, the high-energy excitations can be 
viewed as mobile impurities interacting strongly with the 
underlying TL liquid, and their spectrum can bear sim- 
ilar features as those found in TLL, namely the absence 
of the quasiparticle peak which is replaced by an edge 
singularit y 22 ! 23 : 24 ! 25 It was showr* 2 ^ that the spectrum 
of high-energy excitations with S z = 0, 1 contains an 
edge singularity with a nontrivial field dependence of 
the edge frequency which in the idealized model with 
no interaction except the hardcore constraint is given 
by ujq = (1 — S Z )H . The spectral function of a mobile 
impurity in the TL liquid has been extensively studied 
theoretically ) 26 ! 2 however, to our knowledge, no numer- 
ical results are available to compare with the theoretical 
predictions. An alternative bosonization description in- 
cluding high-energy modes has been proposed recently. 28 
All of this motivates further study of the high-energy 
spectra of Id gapped spin systems in strong field. 

In the present paper, we study the S = 1 Heisen- 
berg chain, which is a paradigmatic example of a Id 
gapped antiferromagnet. There is a large body of nu- 
merical and theoretical results concerning its low-energy 
behavior in strong field o 12 ' 13 : 14 ' 29 which can be used for 
consistency checks. Using the dynamical density matrix 
renormalization group (DMRG) technique, we will show 
that the high-energy part of the spectrum of a S = 1 
chain, placed in a strong external magnetic field H ex- 
ceeding the Haldane gap A, contains edge singularities, 
similar to the well-known infrared singularities in the low- 
energy spectral response. It will be demonstrated that in 
the frequency range uj > A the longitudinal dynamical 
structure factor is dominated by the continuum with a 
power-law edge singularity S^(q — n, uj) cx (uj — ujq) ~ q . 
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The edge exponent a' is found to decrease as a func- 
tion of magnetization, and the edge ujq is shown to 
follow approximately the linear law loq ~ H as found 
earlier^* 2 ! 



II. THEORETICAL PRELIMINARIES 

Consider the 5=1 Heisenberg chain in an applied 
field described by the Hamiltonian 

n = jJ2^n-s n+1 -Hj2s z n , (i) 

n n 

where S„ is the spin-1 operator at site n, the exchange 
constant J will be set to unity and the external mag- 
netic field H is applied along the z axis. In absence of 
the applied field, the ground state of the model is a sin- 
glet, and the lowest excitations are the triplet of Haldane 
magnons separated from the ground state by the gap 
A ~ 0.41 at the wave vector q = tt. The gap closes for 
H > H c = A, and the system acquires a finite density of 
S z = +1 magnons, which can be viewed as bosonic par- 
ticles satisfying the hardcore constraint. It is convenient 
to redefine the momentum so that the minimum of the 
magnon dispersion will correspond to zero. Interacting 
S z = 1 magnons form a TL liquid with the Hamiltonian 
given by 

Ho = y j 'dx{^(d x <p) 2 +K(d x 6) 2 Y (2) 

where ip and 8 is a pair of bosonic fields satisfying the 
commutation relations [<p(x), 0(x')] — iQ{x' — x) (here 
Q(x) is the Heaviside function) and connected by the 
duality relation dtp = vpd x 9. The physical properties of 
the TL liquid are completely characterized by its Fermi 
velocity vp and the TLL parameter K. Both vp and 
K are functions of H, and their behavior for the Hal- 
dane chain has been studied analytically 12,29 as well as 
numerically 13 ' 14 An important feature of the effective 
TLL description of the Haldane chain is the fact that 
the parameter K is always larger than unity (K — > 1 as 
H^H C ). 

The local spin density p(x) (i.e., the density of the sea 
particles) can be expressed as 

p = S z — m + -^=d x p+ const X sin {2kpx+\/^Kw\ , (3) 

where m = m(H) = kp/n is the equilibrium magnetiza- 
tion density at the given field. The operator a* ~ e l V^ s 
creates a kink in the p field shifting it by t/tt and thus 
corresponds to the creation of a S 2 = 1 particle, so 
that at low energies one can establish the correspondence 
S + ~ a) . The low-energy contribution to the transver- 
sal dynamical structure factor (DSF) is thus given by 
S+-{q = tt,w) cx jdxj dte luJt (a^(x,t)a(0,0)), and, since 

(a + (x,t)a(0,0)) oc (v F t 2 - x 2 )^' 2 with r] = —, 



one readily obtains the infrared singularity 

S + ~(q = tt,w) oc a = 2-77 (4) 

with the edge exponent a determined solely by the TL 
liquid parameter K. 

To describe the excitation of "impurity" particles cor- 
responding to the higher-energy magnon branch with 
S z = 0, one can introduce another bosonic field b(x) 
described by the Hamiltonian 

H b = f dx{u tfb+^(d x tf)(d x b)}. (5) 

The scaling dimension of the field b is dim [b] = |. The 
TLL action is Lorentz-invariant and thus dictates the dy- 
namical exponent z = 1, from which one can see that the 
scaling dimension of the mass M is equal to +1, i.e., M 
flows to infinity so that one can neglect the "impurity" 
dispersion described by the second term in ; this argu- 
ment is essentially identical to the dynamical localization 
of a hole moving in an antiferromagnet as discussed in 
Ref. HO- 

We will be interested in the contribution of the excited 
"impurity" particles (^-particles) to the dynamical struc- 
ture factor S zz {q = n,u>), which will be determined by 
^-particles with zero momentum. The local interaction 
between the "impurity" particles and the sea particles (a- 
particles) is proportional to the product of corresponding 
densities b'bp.— This yields, apart from renormalizing 
the energy ujq 1— ► ujq + const • m, two interaction terms: 
one corresponds to forward scattering and is given by 

nL = Uf J dx^=(d x ppb, (6) 

and the other, which is proportional to tfbs'm {2kpx + 
V47r^}, describes backscattering. The backscattering 
term, however, is irrelevant since its scaling dimension is 
1+K, and for the Haldane chain one always has K > 1. 
With the remaining interaction ([S]), the Hamiltonian can 
be diagonalized by means of a unitary transformation 32 

^ —iS / dx6(x)bHx)b(x) KTIt 

U = e J , with (5= L (7) 

Vp^TT 

which changes the phase of b operator according to 

b^b = UHU = be~ m . (8) 

Apart from removing the interaction ([6]), the unitary 
transformation jJJ just renormalizes the threshold en- 
ergy wo, which thus can in principle depend both on the 
interaction and on the applied field. 

Exciting a 6-particle requires "disturbing" the sea of 
a-particles, because the hardcore constraint leads to a 
change in the allowed momentum values when the total 
number of particles is changed^ 3 - To take the hardcore 
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nature of the bosons into account, we can postulate that 
the physical creation operator of the 6-boson can be rep- 
resented as 

^t =6 t a t = bt e i(VW-5)^ (g) 

With that definition, b is the "color-changing" operator 
for an existing sea particle, and in absence of any other 
interaction except the hardcore constraint the excitation 
of a fe-particle would be equivalent to adding another sea 
boson and the only difference would be the additional 
energy cost ujq. 

We are interested in the correlator (i^l(x, t)ipb(®, 0)) 
because it determines the "impurity" contribution to the 
longitudinal DSF 

S**(ir,u)<x J dx J ^e^ t (V' fo t (x,t)Vb(0,0)) (10) 

Since the b quasiparticles are free after the transforma- 
tion, the correlator factorizes. One has (J) [x, t)b(0, 0)) cx 
5(x) e~ luJ<>t , where the Dirac delta-function 5(x) is the 
consequence of irrelevance of the dispersion term in ((5|) , 
so that 

(4(x,t)M0M^e- iuot S(x)/^', (11) 

where 

is the "impurity exponent" . The longitudinal dynamical 
structure factor thus contains an edge singularity 

S zz (n,uj) oc ] —— with a' = l-T}'. (13) 

(w - UJ ) a 

The edge exponent a', in contrast to the low-energy edge 
exponent a, depends not only on the TL liquid parameter 
K, but also on the Fermi velocity vf and on the "host- 
impurity" interaction Uf . In the noninteracting hardcore 
case, which corresponds to Uf = and K = 1, this ex- 
ponent takes the value a! = h, in agreement with earlier 
studies^^^ 

Comparing (TIB"]) and (g]), one can see that the high- 
energy DSF exponent a' = 1 — rf should be consider- 
ably smaller than the low-energy one a = 2 — r\. In 
our approach, this is a direct consequence of the dy- 
namical localization of the "impurity" , which in turn fol- 
lowed from the non-relativistic form of the Tib Hamil- 
tonian ([5]). It is worthwhile to note that the same 
answer is obtained assuming that Tib has a Lorentz- 
invariant form and describes massive field with mass 
ujq and characteristic limiting velocity c. In this lat- 
ter case (P (x,t)b (0,0)) cx Ko(uj r/c), where K a is the 
Bessel function, and r — \J x 2 — c 2 t 2 , and one can show 24 
that the corresponding Fourier transform yields the same 
asymptotic behavior (fT3"|) . 



III. TECHNIQUE 

We use the density matrix renormalization group 
(DMRG) metho d 33 i 34 to calculate the spectral function 
defined earlier. The DMRG works by truncating the 
Hilbert space of the system, based on selecting the opti- 
mal states from the Schmidt decomposition of the lattice 
as split into two (left and right) blocks. The method 
we use is a matrix product generalization of the "correc- 
tion vector method"—, in which we calculate the Green's 
function 

G(q,cj + iS) = (0\Al- 1 - B A q \0) (14) 

E + uj + iS - H H 

at a given frequency u> with Eq being the ground state 
energy of the system and d a numerical broadening, via 
the so-called "correction vector" , 

\c(w + iS)) = (E + w + iS- H)- 1 ^), (15) 

which is determined by a standard linear solver (we used 
GMRES, but biconjugate gradient is also popular). In 
contrast to the traditional dynamical DMRG approach 
where the groundstate |0), Lanczos vector A q \0) and cor- 
rection vector are all determined simultaneously, the ma- 
trix product formulation allows the calculation to be split 
up, improving both the speed and accuracy (for a given 
size of the truncated Hilbert space) of the calculation. 
We first calculate the matrix product approximation to 
the groundstate |0) by a standard DMRG calculation, so 
that we can apply A q exactly to give the Lanczos vector, 
which is independent of frequency. This is used as input 
to the correction vector solver which is trivially paralleliz- 
able over different frequencies, but we can also e.g., use 
a previously determined correction vector from a nearby 
frequency as the initial "guess" vector, additionally ac- 
celerating the calculation. The matrix product formula- 
tion provides good control over the errors, as the residual 
norm || A q \0) — (E + w + iS - H)\c(w + i5)) \\ can be 
calculated exactly. This is prohibitively difficult in the 
standard DMRG approach as it requires calculating the 
matrix elements of H 2 and the usual dynamical DMRG 
approximation H 2 ~ (PHP^) 2 , where P is the projec- 
tor onto the truncated Hilbert space, is inadequate and 
produces an estimate for the residual norm that is many 
orders of magnitude too small. Details of this calculation 
will be provided elsewhere (see also reference [36| for a 
similar approach to this calculation). Given the correc- 
tion vector, the Green's function can be obtained via 

G(q,uj+i5) = (Al\c(uj + iS)), (16) 

from which the spectral function 

S(q,u) = -- lim lmG{q,uj) (17) 

7T (5^0+ 

can be obtained. 
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FIG. 1: Ground state magnetization per site versus magnetic 
field for different system sizes. 



IV. RESULTS AND DISCUSSION 

We start with the low-energy properties of the 
transversal DSF, which are better understood analyti- 
cally, to show the validity of the method. First of all, we 
have to determine the total spin of the ground state de- 
pending on the applied field strength. In Fig. [1] we show 
the ground state magnetization M as a function of the 
magnetic field H; the results for different system sizes 
are consistent, indicating that finite size effects in M(H) 
are small for L > 100. 

The exponent r\ can be obtained from the static trans- 
verse spin-spin correlation function (S x (x)S x (0)) ~ |a;| _,? 
which is easily accessible in DMRG calculations. From 
this the TLL parameter can be calculated via K — 
1/(27/). The values for 77 are shown In Fig. [2] we show our 
values for rj and compare them with analytical results 
based on the effective description in the nonlinear sigma- 
model framework^ as well as with the earlier DMRG 
results. 13 From this r\ we get a first estimate for the edge 
exponent a via the second equation in Eq. Q . 

A typical scan of the low energy-continuum of 
(k = ir, <jj) is shown in Fig. [3] for a 200-site chain. 
The magnetic field applied is H = 0.54 in units of the 
coupling constant leading to a ground state magnetiza- 
tion per site of M/L = 0.05. 

The spectrum reveals a set of peaks at discrete fre- 
quency values, as one expects for a finite-size system. 
The peaks are nearly equidistant (in fact, the distance 
between the peaks decreases slightly with the frequency, 
which is naturally explained by the effective decrease of 
the Fermi velocity as the Fermi sea gets emptied by more 
and more particles being excited). From the energy dif- 
ference of the first two peaks we extract the Fermi veloc- 
ity Vf] its dependence on the magnetic field is in good 
agreement with the results of Konik and Fendlej*^ (see 
Fig. 2]) if one sets their parameter vq to 2.5. This pa- 
rameter has the meaning of a bare spin velocity in the 
nonlinear sigma model description of Konik and Fendley, 
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FIG. 2: Exponent r\ of the static correlation function. Solid 
line: results obtained from analytic calculations by Konik and 
Fendley^, based on the effective nonlinear sigma-model de- 
scription. Open circles: DMRG results for short chain lengths 
L < 80 obtained by Campos Venuti et a/.— Filled squares: 
our results obtained from DMRG calculations using chains of 
length L = 200. 
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FIG. 3: Spectral function S (k = tt,uj) for a 200-site chain 
at magnetic field H = 0.54 leading to a ground state magne- 
tization per site of M/L = 0.05. Here, <5 = 0.01 and m — 300 
states were kept in the calculation. 



and the value vo ~ 2.5 agrees well with the known spin 
wave velocity v « 2.46 of S — 1 Haldane chain^i 

One might be tempted to extract the edge exponent a 
from the dependence of the peak heights in the spectral 
function (Fig. [3]) on the frequency u>. However, this is 
not a very good idea: although all these peaks are nu- 
merically broadened with the same constant S = 0.01 
for a L = 200 site chain, they do not all show the same 
width as it would have been expected. With inceasing 
frequency they get broader. This can be understood in 
the following way: The higher the excitation energy is, 
the more combinations exist to reach it by exciting sev- 
eral particle-hole pairs. In a system with an ideal linear 
dispersion the energies of all those combinations would 
be exactly the same and equal 2?tv pn/L with n being 
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FIG. 4: Fermi velocity extracted from (k = n,u) as a 
function of magnetic field (open circles) . The solid line corre- 
sponds to Eq.(40) of Konik and Fendley^ with the bare spin 
wave velocity vo = 2.5. 

some integer. In a real system the dispersion is not ex- 
actly linear, and those energies become slightly scattered 
around the 2irvFn/L, contributing an additional broad- 
ening to peaks with higher energy. Hence the height of 
the higher energy peaks is underestimated and thus the 
decay is overestimated leading to a value for the exponent 
a being too large. 

There is a much better way to determine the edge ex- 
ponent, namely from the size dependence of the peak 
strength. The height of the first peak at w\ is given by 

S + -{k = Tr,uj 1 ) ocl/wf. 

With u>i cx vp/L we find 

S+-(fc = 7r,wi) cx L a 

for fixed magnetization per site. This is shown in a log- 
log plot in Fig. [5l From calculations for different chain 
lengths L = 50, 100, 150, 200 we find values for a much 
closer to the expected value. 

The values obtained from this L-fit and the expecta- 
tions from a — 2 — r\ are compared in Fig. [SJ The agree- 
ment of both values is good for all M/L (deviation < 5%) 
and gets better for larger magnetizations. Having found 
a reliable way to extract a we turn to the high-energy 
continuum in S zz {k = ir,u)). 

Typical spectra for different magnetizations are shown 
in Fig. [7J The edge frequency shifts to higher energies 
with increasing magnetization. Due to the large en- 
ergy the convergence of the algorithm is worse than for 
the previous continuum resulting in a significantly larger 
DMRG basis. Apart from the expected series of peaks 
these spectra show small additional peaks between each 
two peaks. These are finite-size artifacts that vanish for 
sufficiently large system sizes. 

The Fermi velocity extracted from the upper contin- 
uum spectra is very similar to the one obtained from the 
low-energy continuum as shown in Fig. [5] 



FIG. 5: Peak height of the first peak in S~* (k — ir, to) as func- 
tion of the system size L for different magnetizations. Fitting 
this log-log plot by a straight line yields the edge exponent a. 
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FIG. 6: Edge exponent a obtained from fitting the height of 
the first peaks in the spectral function versus the chain length 
by a power law (full circles). The open squares were obtained 
using Eq|4] 
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FIG. 7: S zz (k = 7T, ui) for different magnetizations. Here the 
numerical broadening is 8 — 0.02 and m = 500 states were 
kept in the reduced DMRG basis. 



6 



2.0 - 



1.5 - 



1.0 - 



B-H upper band 
• lower band 



0.6 



0.8 
H 



1.20 



1.00 



3 0.80 



0.60 




FIG. 8: Fermi velocities extracted from the low-energy spec- 
trum (filled circles) and from the high energy one (open 
squares) . 



FIG. 10: Position of the first peak as function of the magnetic 
field (circles). The dashed line denotes ajn = H as expected 
for the idealized hardcore boson model . 22 ' 23 
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FIG. 9: The edge exponent of the high-energy continuum a' 
obtained from the L-fit of the first peak height for different 
magnetizations. 



FIG. If: Circles: S zz (k — tt,ui) obtained with m — 500 states 
in the DMRG basis. Solid line: Fit by a sum of Lorentzians. 
For fitting parameters see text. 



Like in the low-energy continuum, the peaks get 
broader with increasing energy, making direct power-law 
fit of peak intensity vs frequency difficult. Thus, for de- 
termining the edge exponent a' the same procedure based 
on the size dependence is used as in the previous case. 
From those fits we find the exponent as shown in Fig. [5] 

As a furher check for the consistency of our data, we 
show a fit of the high energy continuum for a magne- 
tization per site M/L = 0.14 calculated in a L = 100 
site chain. We have fitted the spectrum to a sum of 
Lorentzians with different widths and with a power-law 
decay of the peak intensity with the exponent fixed at the 
value a' — 0.79 obtained from the analysis of the size de- 
pendence of the peak heights. The peak widths and the 
true edge frequency ujq were used as fitting parameters. 
The result is shown in Fig.[TTJ Since the small shoulder at 
w « 1.13 is a finite size effect it has not been fitted. The 
values obtained for the peak widths 7, are 71 = 0.0207, 
72 = 0.0236, 73 = 0.0442, and 74 = 0.0820. These widths 



and the independently obtained a' provide an excellent 
fit. The true edge frequency luo — 0.929 is slightly smaller 
than the position of the first peak uii — 0.952. To check 
the validity of this particular value we perform an L — > 00 
extrapolation of tOi(L) which is supposed to yield the 
edge frequency. We find u)\(L — ► 00) = 0.935 which is 
slightly larger than the one found from the fitting proce- 
dure. 

Notably, the exponent a! decreases with increasing 
magnetization unlike the edge exponent of the low-energy 
continuum. With the knowledge of the edge exponent a', 
the Fermi velocity Vf and the static exponent i] we can 
extract the effective '"host-impurity" interaction Uf us- 
ing Eqs. (|12|) . (fl3| . The values of Uf extracted in that 
way lie between 7.7 and 8.3, which can be deemed approx- 
imately independent of the magnetization (within ±5% 
accuracy) . 

The position of the first peak toi (which serves as a 
crude approximation for the true edge frequency ujq) as 
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a function of the magnetic field shows a clear linear de- 
pendence (see Fig. [TU)) . It should be noted that a linear 
dependence u>o = H has been obtained in the idealized 
hardcore bosonic model ) 22 ! 23 as well as in a bosonisation 
calculation 2 ^ for a S — \ ladder. The observed depen- 
dence is quite close to loq = H as seen in Fig. |TD] 

V. SUMMARY 

We have shown that the high-energy spectrum of a 
S = 1 Heisenberg chain in a strong external magnetic 
field H exceeding the Haldane gap A contains edge sin- 
gularities, similar to those known to exist in the low- 
energy spectral response. It is found that in the fre- 
quency range to > A the longitudinal (with respect to 
the applied field) dynamical structure factor is dominated 
by the power-law singularity 5" (q = n, lu) oc (lu — LUo)~ a ■ 
It is shown that the edge exponent of the high-energy 
continuum a' decreases with the magnetic field, consis- 



tent with theoretical expectations. The edge frequency 

luq is found to increase linearly with the magnetic field, 

uiq ~ cH, the coefficient c being close to 1, which agrees 

i — ^ii — ' — ir — i 

well with the predictions of Rcfs. 22.23 24. The existence 
of power-law continua with edge singularities at high en- 
ergies should be a quite general feature, common for all 
one-dimensional gapped spin systems in high magnetic 
fields and directly related to the Tomonaga-Luttinger liq- 
uid character of the ground state at H > A. 
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